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COMPUTATION OF VISCOUS BLAST WAVE FLOWFIELDS 

Christopher A. Atwood 

The study of the effects of blast-wave impingment upon 
vehicles and structures is of practical consideration in the 
determination of their survivability. Unfortunately, the 
experimental study of the blast-wave/target interaction problem is 
expensive, and detailed flowfield quantities are difficult to obtain. 
However, recent advances in computational fluid dynamics may 
provide a means, complementary to experimental studies, by which 
the timely design of effective configurations can be found. 

A method to determine unsteady solutions of the Navier- 
Stokes equations has been developed and applied. The structured 
finite-volume, approximately factored implicit scheme uses Newton 
subiterations to obtain the spatially and temporally second-order 
accurate time history of the interaction of bast-waves with 
stationary targets. The inviscid flux is evaluated using 
MacCormack's modified Steger-Warming flux or Roe flux difference 
splittings with total variation diminishing limiters, while the 
viscous flux is computed using central differences. The use of 
implicit boundary conditions in conjunction with a telescoping in 
time and space method permitted solutions to this strongly unsteady 
class of problems. Comparisons of numerical, analytical, and 
experimental results have been made in two and three dimensions. 
These comparisons revealed accurate wave speed resolution with 
nonoscillatory discontinuity capturing. 

The simulation of the inviscid blast-wave problem has been 
studied in the past by Champney, Chaussee, and Kutler; by Mark and 
Kutler; and by Lohner and Yee. The addition of viscous effects has 
been studied in two dimensions by Bennett, Abbett, and Wolf; by 
Molvik; and by Hisley and Molvik. The purpose of this effort was to 
address the three-dimensional, viscous blast-wave problem. 

Test cases were undertaken to reveal these methods' 
weaknesses in three regimes: 1) viscous-dominated flow; 2) 

complex unsteady flow; and 3) three-dimensional flow. Comparisons 
of these computations to analytic and experimental results provided 
initial validation of the resultant code. Additional details on the 
numerical method and on the validation can be found in Appendix A. 
Presently, the code which has been released is capable of single 


zone computations with selection of any permutation of solid wall 
or flow-through boundaries. 

The disadvantage of these characteristic-based schemes is 
their expense: 86 ps/cell/iteration at 140 MFLOPS on a single head 

of the CCF ray Y-MP/832. However, memory use is a conventional 38 
words/cell. Incorporation of a zonal methodology, a turbulence 
model and time metrics would increase the versatility of the 
present code. 
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Abstract 

A numerical method to determine unsteady solu- 
tions of the laminar, perfect gas Navier-Stokes equa- 
tions has been developed. The structured finite- 
volume, approximately factored implicit scheme uses 
Newton subiterations to obtain the spatially and tem- 
porally second-order accurate time history of the in- 
teraction of blast- waves with stationary targets. The 
inviscid flux is evaluated using either of two upwind 
techniques, while the viscous terms are computed by 
central differencing. Comparisons of numerical, ana- 
lytical, and experimental results are made in two and 
three dimensions. The results show accurate wave 
speed resolution and nonoscillatory discontinuity cap- 
turing. 

Nomenclature 

A, B } C inviscid flux Jacobians 

c speed of sound 

c T specific heat at constant pressure 

c* specific heat at constant volume 

e total energy per unit volume 

E , F , G flux vectors 

F flux tensor of second order 

h enthalpy per unit mass 

i,j,k Cartesian unit vectors 

J coordinate transformation Jacobian 

Jb coefficient of thermal conductivity 

M Mach number or viscous flux Jacobian 

p static pressure 

Pr Prandtl number 

q velocity magnitude 

q heat transfer vector 

Q vector of dependent variables 

R specific gas constant 

Re Reynolds number 

r position vector in physical space 

3 entropy function 

s outward-directed surface normal 

t time 

T absolute temperature 

T“ l ,T matrices of left and right eigenvectors 

u, v , w Cartesian velocity components 
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u,v,w 

contravariant velocities 

V 

volume 

x,y,* 

Cartesian physical space coordinates 

0 

compression parameter 


ratio of specific heats 

€ 

internal energy per unit mass 

c 

bulk coefficient of viscosity 

K 

ratio of coefficient of thermal conductivity 
to the specific heat at constant volume 

A 

second coefficient of viscosity 

A 

diagonal matrix of eigenvalues, A, 

P 

dynamic or first coefficient of viscosity 

(,V,C 

curvilinear space coordinates 

P 

density 

T 

computational temporal coordinate 

T ij 

viscous stress tensor 

<f> 

flux influence parameter 

Superscripts 

i 

wave family 

m 

subiteration level 

n 

time level 

Subscripts 


rji C direction indices 

NS 

Navier-Stokes 

T 

total or stagnation quantity 

x,y,z 

partial with respect to Cartesian 
coordinate 


partial w.r.t. curvilinear coordinate 


Introduction 

The study of the effects of blast-wave impingment 
upon vehicles and structures is of practical consid- 
eration in the determination of their survivability. 
The experimental study of the bleist- wave/target in- 
teraction problem is expensive, and detailed flowfield 
quantities are difficult to obtain. Recent advances in 
computational fluid dynamics may provide a means, 
complementary to experimental studies, by which the 
timely design of effective configurations can be found. 

The simulation of the inviscid blast-wave problem 
has been studied in the past by Champney, Chaussee, 
and Kutler , 1 Mark and Kutler , 2 Lohner , 3 and Yee . 4 
The addition of viscous effects has been studied in two 
dimensions by Bennett, Abbett, and Wolf , 5 Molvik , 6 

t 
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and Hisley and Molvik . 7 The purpose of this effort is 
to address the three-dimensional, viscous blast- wave 
problem. 

The techniques developed here utilise total varia- 
tion diminishing (TVD) upwind and upwind-biased 
schemes to resolve the discontinuous flow features 
without the oscillations prevalent in the more con- 
ventional of the central difference methods. Wave 
speeds are resolved adequately at large Courant num- 
bers through the use of time conservative differencing 
and Newton subiterations to retain accuracy. 

The following sections provide background material 
on the assumptions made, the analysis of the implicit 
scheme, and the flux evaluation methods. The exten- 
sions of previous work are also noted. Finally, the 
methods are applied to five test cases, including com- 
parisons of numerical, analytical, and experimental 
results. 


The Governing Equations 

The Navier- Stokes Equations 

The Navier-Stokes equations may be expressed in 
integral conservation law form, coupled with the con- 
tinuity and energy equations as 

(1) 

where body forces have been neglected and the cell 
volumes are time invariant. Here V is the volume 
of an arbitrary fluid packet, f is the flux tensor of 
second order, and ds is an outward directed normal of 
a differential surface area. The vectors may be written 
in Cartesian coordinates as 


Q = \p,pv,pv,p w,e] T 


Ens = 


pu 

pit 1 +p+ r„ 
puv + r,, 
puw + r„ 

(e + p)tt + r„ii + T,,t> + w + q. , 


Fits = 


pv 

pVU + Ty t 

pv 7 +P + T„ 

pvw + Tyt 

(e + p)v + TytU + TyyV + Ty, W + 


Gus = 


pw 

pU)U + T sa 

fnvv + r,y 
pw 7 +P+T„ 

(e + p)w + r lx u + r^v + r„ w + q t 


where each flux can be partitioned into inviscid and 
viscous portions. The viscous portion is composed of 


terms: 
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The total energy per unit volume is related to the 
internal energy per unit mass by e = pe + pq 7 / 2. The 
perfect gas equation of state, p =r pRT % completes 
the system. In addition, for thermally and calorically 
perfect gases the specific heats are constant, leaving 
e = c,T, and h = CpT. The ratio of specific heats is 
taken as 7 = = 1.4, and the thermodynamic states 

are related using 

*r = fc±^, + + »’ + «■) 

Fourier’s law for heat transfer by conduction is as- 
sumed; hence, the heat transfer can be expressed as 


q = -*VT = -(?,i + 9 „j +?.k) 

(9e 8e 8(\ 

V^ l+ V ^ J 


= —/c 


where k = Jk/c, = 7 pfPr. The Prandtl number for 
air is fixed at Pr = 0.72. 

The relationship between the first (/ 1 ), second (A), 
and bulk (£) viscosity coefficients is C = §M + A. The 
bulk viscosity coefficient is set to sero in accordance 
with Stokes’ hypothesis, resulting in A = — 2/3/x. Vis- 
cosity is related to the thermodynamic state using 
Sutherland’s formula: 

C t Ti 

T+C, 

where C\ and Cj are specific to the gas in question. 

A solution to this system of nonlinear partial differ- 
ential equations may be obtained numerically over a 
discretised domain, which in turn is frequently trans- 
formed into a convenient computational space. 


Transformation to Curvilinear Coordinates 

In order to adequately resolve the solid bound- 
ary/fluid interaction, it is common to transform the 
above equations into curvilinear coordinates which 
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can be body-conformal. Specifically, the body is con- 
strained to lie at a constant £,17, or £ value. For a 
stationary grid, this transformation can be expressed 
as 


T = t, £ = £(x t y, z), 7 ) ~ tj(x, y, z), ( = <(*, y, z) 
Application of the chain rule of differentiation yields 

a a a a 

3 * “ ** AC + * dv + C * 

with similar expressions for the partials with respect 
to y and z. The inverse transformation gives 

a a a a 

Again, expressions can be found for the tj and £ par- 
tials in a like manner. Represented in matrix form: 

£* W* C* 

ty Vy Cf 

6 v* (* J 

T 

and for the inverse transformation, 

*( V( H " 

Vi *n 

*< V< z < . 

-v 
T- 1 

Combining the use of T = (T _l )~ l and finite volume 
metrics, such as those described by Vinokur, 8 leads 
to a scheme which is freestream-preserving because of 
the telescoping property. Hence, if the surface nor- 
mals to a constant £ , 17, or £ plane are defined respec- 
tively as 

- V»'+J i+ + 

= ^(rr - r 4 ) x (r 8 - r s ) 

*>+i = *«J+J i + Vi+jj + , *J+J k 

= i(r 7 -r 2 ) x (r s -r e ) 

S *+J = *«,*+|» + *y,k+jj + V*+$ k 

= ^(r 6 -r 8 )x(r 6 -r 7 ) 

where the index convention is shown in Fig. 1 , then 
the metrics can be formed as 

<• = J hhH~Vi z n) = 

— X n Z() = y*y,«+J 


8 1 
as - 
a 

3 y 

a 

Tz J 


r a 1 

a? 

a 

a 

K J 


a i 
a? 
a 

3 y 

a 

a? J 



6 = J(*.»yc -*<*>) = y^.i+i 

ij, = J(y<z ( -y(Z ( ) = 

»j, = J(*<z< -*<**) = 

Vi = -*{»() = ^'ij+l 

C* = J(y(*v-yn*() = 

C, = J(*,z { - z ( z,) = yVk+| 

<* = J(*ey,-*i,»E) = 

These metrics represent the projections of the cell face 
normal into (x,y, z) space. The faces of the hexahe- 
dron exactly enclose the discrete control volume, i.e., 
no gaps are permitted at the edges. Finally, the Ja- 
cobian of the coordinate transformation is equivalent 
to the inverse of the volume, as related by 

£ _ a(z, y, *) 

j d(C,V,C) 

= *{(y->*< - y<*n) - *ij(y« z < - y< z () 

+*<(yf*i - Vv z () 

= y = +»*_$) -(ry-ri) 

Utilizing these metrics in the application of the 
chain rule to Eq. (1) and subsequent simplification 
yields 

Q' r + E' ( + F; + G't = 0 

where 

Q' = QV 

E'nS = (Essix + Fssiy + G N si,)V 

— [EnS*z + FnS» y + G}ts*i |< 

Efts = (EnsVx + EitsVy + GffSVt)y 
= [Ens»z + Fn S»y + Gns*i I j 
G'ns = {EssCx + EitsCy + GffsCt ) V 

= [-Ejyrs** + -PWsZy + Gj\r5S x |fc 

Separating the inviscid and viscous portions of the 
flux vectors, in the £ direction E' NS = E' + E' v , then 

’ pU ' 

fmU + ixP 

E' = V pvV + i y p 
pwU + isP 
. (e+P)V . 

Here the contravariant velocity is U = u£ z +v£ y +w£ t , 
without metric normalization. The viscous flux can be 
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represented as 


The remaining variables follow: 


E' 9 =V 


0 


+ Tyzty “I" T *zts 
+ T yyty Tiy^i 
r 9M (z + v£y + Txsil 

(«e' a -f ve* 3 + we't) + {q*tx + qyiy + ^ 6 ) . 


where the viscous stress terms axe evaluated by again 
invoking the chain rule) and the flux in the rj and £ 
directions are found similarly. The results presented 
herein are implemented using either the thin-layer or 
the full viscous term treatment. 

The widespread use of the thin-layer approxima- 
tion, first implemented by Steger, 9 can be justi- 
fied from either physical or algorithmical arguments. 
Physically, the neglect of all diffusion processes par- 
allel to the body is similar to that used in boundary- 
layer theory, albeit not as restrictive. Hence, when the 
viscous effects are confined to thin regions which fall 
along a constant £ , 77, or £ plane then this assumption 
is valid. Algorithmically, the banded matrix struc- 
ture used in multidimensional algorithms which se- 
quentially solve a set of unidirectional problems can 
include only these thin-layer terms implicitly. This 
thin-layer flux in the rj direction, assumed to be the 
body normal coordinate, is expressed as: 

0 

miUjj + m A v v + msw v 
m+Uij + m 7 v v -f me w v 

7715 U, + + TTlsWr) 

miliUtj + m 7 vv v + m&ww v 

+ m 5 (tin;,, + wu,) 

-f- wv v ) -I- K€r,{T]l + + r£) m 



where the (“) denotes an arithmetic mean value and 


TOx 

= p(jvl + vl+vlj 

| . ™4 

m 2 

= #*fn2 + !«£+»£) 

| . 

7713 

= p (vl + vl + j 

| . 


P 

3 ** 
P 

jW* 

P 


These viscous flux terms may be found for the remain- 
ing spatial coordinates as well. The results presented 
here are implemented using either the thin-layer or 
the full viscous term treatment, as required by the 
the flow physics. 


P = P/Pr.f, P = P/Pr.f, e = e/(Pr./t*?«/) 

U = i/dr.fj V = v/Or'J, W = w/Ur t) 

t = tUr.f/L, T = p/p, H = A. /At./ 

The Reynolds number resulting from this procedure 
is Re ~ Pr t f Lurtf/jlrcf , where the (“) denotes a di- 
mensional quantity, and the (<») denotes the quies- 
cent conditions surrounding the target before primary 
shock arrival. 


The Numerical Technique 


The development of the scheme will be described by 
discussion of the first-order terms, following which the 
higher-order extensions will be outlined. The scheme 
as expressed for a cell which has a mean flux value on 
each of the six sides is 



+ 

+ 





+ / , / . ( G U*+i = 0 


In discrete form, after dropping the primes for conve- 
nience, the governing equations can be written as 




+ <«r*i+* - ® 


where n denotes the time level in this implicit repre- 
sentation, and A£, At;, and AC are set to unity for 
convenience. 

These flux terms may be evaluated using a tech- 
nique which may be broadly classed as either central 
or upwind. The latter technique is chosen for this 
study for the desirable numerical properties, such as 
diagonal dominance of the flux Jacobian, and for the 
physical dependence on zones of influence which are 
inherent in upwind schemes. 


Nondimen sionalizat ion 

The governing equations may be nondimensional- 
ized by the choice of a length scale, denoted and 
reference values of p, v, and p such as 

Prtf — Poo* — VPoo/ Poo* Prtf ~ Poo 


Upwind Schemes 

Upwind schemes bias the derivative evaluations re- 
quired to determine the flux across fluid cells accord- 
ing to the sign of the eigenvalues. In this manner 
these methods bring the physics of the hyperbolic sys- 
tem, the unsteady Euler equations, into the numerical 
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solution process. To facilitate the implementation of 
these upwind schemes, the eigensystem is determined. 
The similarity transformation which diagonalizes the 
unsteady, inviscid, gas-dynamic equations, shown by 
Warming, Beam, and Hyett, 10 is outlined as follows 

E = AQ = TAT~ l Q 

where the rows of T -1 are the eigenvectors and 

A = A + + A" = diag [u, U,U,U + c,U - c] ||s|| 

using normalized contravariant velocities. The eigen- 
values can be split according to, among other split- 
tings, their signs: 

_ A ± |A| 

2 

Two upwind schemes are implemented here to com- 
pare the results which may be obtained with either of 
the techniques. The early portion of this discussion 
will be shown unidimensionally; the multidimensional 
extension will be outlined towards the end of this sec- 
tion. 

Flux Vector Splitting 

The shock-capturing scheme developed by Steger 
and Warming 11 revisited the classical characteris- 
tic procedures. They found that the Euler equa- 
tions possessed the property of homogeneity of de- 
gree one for the equation of state used here, meaning 
E(aQ) = a£(Q). For a vector with this property 
E = AQ and consequently the flux vector can be split 
into two parts, each physically corresponding to the 
right and left moving waves. This technique resulted 
in the flux being represented as a combination of the 
subspaces associated with the positive and negative 
eigenvalues, expressed as 

E = T(A + + A~)T~ l Q = {A+ + A~)Q 
= F + + £T 

where T and T _1 are the right and left eigenvectors, 
respectively. The flux across a cell face can be deter- 
mined by 

E i + i = K+i + E 7+i 

= *t+ 

Because the Jacobian at i + | is dependent on two 
states, this solution method now diverges from the 
original flux vector splitting. The treatment of this 
Jacobian is shown in a following section. Linearization 
in time can be performed for one flux as follows, the 
extension to the remaining flux terms is omitted for 
brevity. At the n + 1 time level, 

K+i = m + ): HQi+ i +{A~)7}iQi+i 

= {A + ): +i 6 Qi + 


where the implicit change in the dependent variables 
is given by SQ = Q n+l — Q n . Note that the Jacobian 
matrices are frozen at time level n. The remaining 
flax, may be obtained similarly. 

To assess the effect on stability of this type of 
linearization, a procedure given by Barth 13 is ap- 
plied to this method. Using the semidiscrete form 
dQi/di = —Rit then 



Using frozen Jacobian matrices, the method can be 
linearized as follows: 



where for the first-order subset the Jacobian is a block 
tridiagonal matrix. The blocks along the I th row are 



This scheme is inherently conservative in space be- 
cause of the telescoping property of the finite-volume 
formulation; analysis of this scheme reveals that it 
is also conservative in time, thus allowing the use of 
large Courant numbers. A demonstration of this anal- 
ysis proceeds by writing the scheme as 

-L + A"J SQ = -R n = - A n Q " 

hence 

Q n+1 + AiA n Q n+1 = Q n (3) 

In order for the scheme to be conservative in time over 
a periodic domain, the global average of a solution 
must remain constant for all time, i.e., J^i=i Q? = 

Qi = ELx < ?" +1 - Hence, when Eq. (3) is 
summed across the domain, the result is that the 
columns of A n must sum to zero. For the scheme out- 
lined here, this can be verified with some effort. Use of 
these linearizations obviates the need for convergence 
at the subiteration level to obtain time conservation. 

Flux Difference Splitting 

Flux difference splitting methods are based on 
the Riemann problem, solved exactly by Godunov in 
1959. 18 The Riemann problem is composed of m -f 1 
piecewise constant states separated by m wave fami- 
lies. The waves include shocks, contact surfaces, and 
rarefaction fans. For each of the Riemann problem 
cells, the transition of the dependent variables is a 
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function of a parameter family. The solution can be 
found once these transition states are known. Approx- 
imate Riemann solvers simplify the numerics of the 
problem by eliminating the iterative process required 
to find the intermediate states. Figure 2 shows a 
schematic of the Riemann problem with the piecewise 
constant states separated by the appropriate wave 
families. The flux through the cell face is 

*<+* = 2 ZI [* + ^ i+1 " (A ^+i “ Af wp] 

where \E\ = |A|Q = ( A + - A~)Q = [T(A+ - 
A^)T- l )Q. The flux differences associated with the 
■f and — traveling waves are 

A E* +i = (TA*T~ l ) i+i (Qi+i - Qi) 

A E‘ +i = (TA~T~ l ) i+ ^(Qi+i — Qi) 

Again utilizing the semidiscrete form dQij&i = — i2», 
then 


where = IT + A n 6Q and <? n+l = Q n + 

( dQ/dt) n At . The flux through the remaining faces 
are determined similarly. This scheme can also be 
shown to obey the criterion for conservation in time. 

Roe Averaging 

In order to determine the Jacobian at the cell face 
* + some function, A = A(QliQr)> must be as- 
sumed, where the subscripts indicate left and right 
states. The location where this flux evaluation occurs 
is among the discrepancies between finite-difference 
and finite- volume schemes. The Jacobian form used 
here is attributable to Roe, 14 which provides an ap- 
proximate solution to the Riemann problem. This Ja- 
cobian is created through the use of a parameter vec- 
tor composed of a geometric-like mean of the states. 
The more obvious arithmetic Jacobian forms, such as 
A = \(Al+Aj 1 ) or A = are not con- 

servative forms. Conservative Jacobian forms satisfy 
A(Ql, Qr)(Ql - Qr ) = El — Er . Stated explicitly, 
the Roe averaging operation is 


Ri = 2 ^ [(^+1 - Ei-i) - (A|S|« + J + A|£| M )] 

This method can be linearized using a procedure simi- 
lar to that described previously in Eq. (2). The blocks 
of the i th row are now expressed as 


dRi _ 1 ( 

dQi _r 2Ax ^ ^ _1 + dQi _ a ) 

dRi _ 1 ( 0A|£| i+i *A| EU.A 

dQi 2A* ^ 8Qi + dQi ) 

dRi _ 1 ( ^A|f?| j+i \ 

dQ i+1 2Ax ^ ’ +x dQ i+1 ) 

Substitution of the flux difference splitting expression 
yields 


dA\ Elj+i _ Q 

dQi+1 dQi+1 


= |A| i+ j + 


[l-^li+ ““ Q«) ] 


dQi+i 


(Qi+i-Qi) 


These true Jacobians are expensive to compute, and 
the simplification to approximate Jacobians is made 

dQ%+\ 

Utilizing these approximate Jacobians, the linearisa- 
tion proceeds as 






^[m?+kv } w< 

+W«-W +1 )*«m]+3V t 


yJpLPR 

V\LyfPL + UiRy/PR 

>fpZ + yfpR 

PL*iL + Pjvil, + ViR) + PRViR 
PL + 2p + PR 

Pl{^t)l + p{{^t)l - I - (^t)r) + pr{^t)r 

PL + 2p + PR 

where a (”) denotes a Roe averaged quantity and the 
latter forms are presented as inexpensive alternative 
expressions. Substitution of (i) for L and (* + 1) for 
R allows the evaluation of the Jacobian at the inter- 
mediary cell face. For the flux vector splitting case 
described earlier, MacCormack 15 has found this aver- 
age helps to alleviate excessive numerical dissipation 
in regions dominated by viscous effects. Roe aver- 
aged values are utilized throughout the development 
presented here. 

Higher-Order Extensions 

Spatially first-order methods frequently provide in- 
adequate resolution of the flowfield. However, the 
methods discussed above can be extended to higher- 
order spatial accuracy by modification of the right 
hand side. In order to assist in the preservation of 
well-behaved solutions near the discontinuities admit- 
ted by the strong conservation law form of the Euler 
equations, a total variation diminishing technique is 
implemented. If the total variation of a solution is de- 
fined as TV(u) = X/fl -00 |uf+i ~ then a solution 
which follows TV{u n + l ) < TV{u n ) is TVD. The TVD 
constraint can be shown to result in diagonal domi- 
nance, allowing the use of relaxation schemes. In this 
manner the scheme may be extended to higher space 
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accuracy thoughout the smoothly varying regions of 
the field, reducing the accuracy in localities of high- 
gradient and extrema in order to obtain sharp and 
oscillation-free resolution. These methods are rigor- 
ously applicable only to scalar nonlinear equations or 
a system of linear equations in one spatial dimension. 
Application of these schemes to multidimensional sys- 
tems of nonlinear equations are generally not TVD. 
Moreover, it is not clear that the higher-order accu- 
racy of the unidimensional problem is retained in mul- 
tidimensional cases. However, the results which can 
be obtained demonstrate the usefulness of the tech- 
nique. 

Of the several methods which fall into the TVD 
domain, 4 the technique implemented here is one at- 
tributable to Chakravarthy and Osher, 16 the develop- 
ment of which follows for completeness. In this formu- 
lation, the higher-order flux can be expressed as a sum 
of a first-order flux, denoted 22,+j, and a flux correc- 
tion term. The flux correction terms are determined 
by first computing the flux differences across the m 
wave families mentioned previously. Subsequent lim- 
iting of these flux differences and summation across 
the wave families results in the higher-order flux. This 
flux is expressed as 

^>+i = E i+i 
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where (~) and (~) indicate a quantity that has been 
limited, j is the index denoting the wave family, and 
i is the index assigned to a cell center. Using the no- 
tation of V for the rows of the left eigenvector matrix, 
T -1 , and r* for the columns of the right eigenvector 
matrix, T, then the measure of the change in the de- 
pendent variables is 

"<+$ = ■ (<?•+! “ Qi) 

The measure of the change in the flux is defined as 

<r>= AV = (A* + +A'-)tt' 

the eigenvalues being split as shown previously. The 
limited counterparts of these values are obtained as: 


= minmod 


<r i+ j = minmod 


j = minmod 


d\_ | = minmod 



This limiter returns the argument of smaller magni- 
tude when the signs are equal, and returns zero when 
the arguments are of opposite sign. This procedure 
effectively adds dissipation locally in regions of high 
flux gradient and at inflection points. In this manner, 
monotonicity is preserved by preventing the creation 
of new extrema while preserving the global accuracy 
of the solution. While formal accuracy estimates are 
difficult to ascertain because of the nonlinear applica- 
tion of limiting to different wave families, numerical 
experiments have demonstrated that the global accu- 
racy of the underlying scheme is preserved. 17 

The compression parameter, /3, is restricted accord- 
ing to 1 < 0 < and the limiting operator is given 
as 

minmod(x,y) = sign(x) (max{0, min [|x|, ysign(x)]}) 

The compression parameter reduces the amount of 
dissipation added, the range being bounded by ac- 
curacy and TVD constraints. Finally, the limited flux 
difference values are expressed as 
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This asymmetric limiter is designed to modify the 
fluxes only in the rapidly varying portions of the 
flow, where nonphysical oscillations are likely to oc- 
cur. Since these high-gradient regions are confined 
to thin regions, the dominant solution domain is dif- 
ferenced in accordance with the underlying scheme. 
Variances in the value of the compression parameter 
allow the fluxes to be limited for different gradient 
levels. This implies that use of 0max will cause the 
limiting action to be taken only in the high-gradient 
regions, and lower values of /3 will result in limiting 
for commensurately lower flux gradients. The variety 
of schemes which can be obtained using this technique 
are shown in Table 1. 


Table 1: Summary of schemes. 


<t> 

Unlimited Scheme 

0max 

2 nd order TE 

-1 

Fully upwind 

2 


i 

3 

Nameless 

5 

2 

|(Ax) 2 / xxx 

0 

Fromm’s 

3 


1 

3 

Z Td Order 

4 

0 

1 

3 

Low TE 2 nd Order 

5 

-£(A *)’/.« 

1 

Central 

oo 



Here TE = (| — ^)(Ax) 2 / xxx /4 defines the leading 
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term of the truncation error for the unlimited form 
of the schemes. Local metrics have been used in the 
above method to maintain reasonable computational 
efficiency, a satisfactory approximation for grids which 
do not contain rapid variations. 


Viscous Terms 


The viscous terms are treated through central dif- 
ferencing about the cell faces. The explicit terms are 
conventionally differenced after chain-rule expansion, 
inclusive of the cross terms if these diffusion processes 
are deemed significant. The left hand side does not 
include these cross terms, and the resultant viscous 
Jacobian, employing V = [p, u, v, tt>, e] T as the primi- 
tive variable vector, is 


M = - 


0 0 0 0 0 

0 mj 3 M* x s y /3 p***j/3 0 

0 77133 T7I33 P*y*x/3 0 

0 77124 77134 77144 0 

0 77152 77I53 77154 771 5 5 


77*22 =M( 3 ** + <y + , 77133 — /*(*» + 2*1 + 

m AA =n(a\+ s\ + , m ss = *{»l +*l +*J) 


77152 = U 77122 + V 77123 + ttff7l24 

771 5 3 = 1177132 + V77133 + 1^77134 

77154 = 7*77142 + V 77143 + 7^77144 


Now, expansion of the block structure gives 




+B ij+ }.* 


| "vT^.+f j.k = 


where N = dV/dQ and only the thin-layer terms in 77 
are shown here. 


Factorisation 


The extension of the unidimensional techniques 
given above is accomplished through dimensional 
splitting. The method used here is that of Yanenko, 18 
where the factors are chosen in the 77, and ( direc- 
tions. Expressing the three-dimensional equations in 
compact notation as 




then the factorisation procedure yields 


(/ + ^)(/ + ^B)(/ + |r C )«< ? =Ao 


This system can be solved sequentially through the 
use of intermediary steps without loss of time accu- 
racy. Although alternating direction implicit schemes 
of this type offer advantages of vectorixation, the sys- 
tem is solved as a sequence of unidimensional prob- 
lems, hence limiting the size of the time step due to 
stability restrictions. 1 ® The use of this technique here 
is justified by the requirement of adequate time his- 
tory flow resolution, thus imposing an additional con- 
straint on the maximum feasible time step. Appli- 
cation of a line Gauss-Seidel method to the second 
test case confirmed this hypothesis. This relaxation 
method offered a slightly increased stability range, but 
not enough to warrant its additional expense. Ad- 
ditionally, since for the factored scheme the flux ex- 
change occurs at the same time level, the technique is 
conservative, even when convergence at the subitera- 
tion level is not attained for each time step. 


Newton Iterative Technique 

Reduction of the linearization and factorization er- 
rors is achieved by a Newton iterative method of the 
type described by Rai and Chakravarthy 30 and Rogers 
and Kwak 31 is utilized, albeit with the addition of al- 
lowance for a varying step size.® Assuming that the 
initial guess lies within the radius of convergence, the 
right hand side is converged to an arbitrary accu- 
racy while holding time fixed. Since the right hand 
side includes the higher-order difference representa- 
tions of the Navier-Stokes equations, these errors are 
eliminated at convergence. The method is discussed 
below where m is the iteration index and n is the 
conventional index denoting time level. Discretizing 
Qt + E* = 0 gives 


^ (gu+l.m+1 _ gn+l,mj 
- g"+*.m+l + J_(gn + l,m _ gnj 
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+ 


_ ^ jQn + l.m + l _ gn,m+ 1 _ (Q n + 1 * m „ Q n 

Defining 6Q 9 = Q n+l * m+I - Q n+l ' m , then 

J£Q' = Ar<2" +l,m+l - (Q n + 1 » m - Q n ) 

= -Ar££ +1 ’ m+1 - (Q n+1 * m - Q n ) 

Linearization at iteration level m - h 1 gives 


to attain rapid convergence. The left hand side allows 
the use of large time steps by relieving the Courant- 
Friedrichs-Lewy stability constraint. Deferring the 
question of uniqueness, if a set of dependent variables 
is found such that the right hand side is satisfied, then 
this field is a solution to the discretized equations re- 
gardless of the approximations made to arrive at that 
set. 




/0££ +1 \ 

\ dQ ) 


6Q' 


where the flux Jacobian has been frozen at iteration 
m. Substitution yields 

I6Q' - -ArE^ +1 ’ m -AT^-A n+hm SQ'-(Q n+l ’ m -Q n ) 


Rearranging results in 

Jj+ ^I J 4 n+1 ' m j 6Q' = -(Q n+1 ' m -Q"+Ar£^ +lim ) 

which reverts to the standard noniterative form when 
no subiterations are taken, as can be seen by substi- 
tution of n for n + 1, m. 

The temporally second-order accurate representa- 
tion is found by extension of the above procedure. 
Using a three-point backward time stencil, 

Qt = C 0 Q n+1 + CrQ n + C 7 Q n ~ l 


where 


Co 


C j 


IzZ c = _ 

(1 - <r)Ar 2 + An ’ 1 (1 - <r)AT 2 -f An 

-1 

(1 — <t) Atj + A T\ 


and <r = (1 + Ati/Atj) 2 . The elapsed time between 
the n — 1 and n time levels is given by Arj. and the 
step size between n and n + 1 is A r 3 . Finally, 


U + 


Cq Ax 


-A n+1 ’ m 


6Q' = 


-(Q n+1,m + ^Q n + %-Q" 

Cq Oq 



which again reduces to the special noniterative case. 
The use of a variable time step allows the solution 
to progress using a constant Courant number, possi- 
bly preventing inadvertent divergences. Higher-order 
accuracy in time may be determined by extension of 
the above technique, albeit with additional memory 
requirements. 

The assertion that this technique reduces the fac- 
torization and linearization errors is substantiated as 
follows. The right hand side of the method contains 
the discretized governing equations in their pure form, 
that is, without the numerical approximations utilized 


Grid Generation 

Grid generation is a significant portion of the ef- 
fort spent in obtaining the flowfield about any rea- 
sonably complex geometry. A structured approach is 
utilized in this study, with the generalized grids gen- 
erated using codes written by Steinbrenner, Chawner, 
and Fouts, 33 and Atwood and Vogel. 33 The cases cho- 
sen here were, to a certain extent, driven by the an- 
ticipated time costliness of the surface definition and 
mesh generation. 

Boundary Conditions 

The block implicit boundary conditions are imple- 
mented in a manner consistent with the previously de- 
scribed flux split linearization described earlier. The 
inviscid and viscous impermeable wall conditions are 
prescribed similarly to those given by MacCormack. 34 
Although the following procedures are presented for a 
cell face which lies along a constant £ plane, the pro- 
cedure may be generalized for application to any cell 
boundary. 

The inviscid, impermeable boundary condition is 
described for a pair of cells between which the surface 
lies, depicted in Fig. 3. In the following discussion, 
the cell above the wall will be denoted by subscript 2, 
the cell below the wall by subscript 1. At cell centroid 
2 the velocity is expressed as 

V* = ui + vj + wk\i 

and at the inviscid wall the surface normal is 
*wall — p(i.»+ *yj "t" *xk)| wa /j 

= SyJ + ?*k)| wall 

hence the velocity component normal to the wall is 
= H^nallS-.^ 

= (us* + vty + + SyJ + s,k) 

Since, = Vt + then the tangential velocity com- 
ponent is 

= [u - + vty + 

+[r - ? y (us x + vt y + u)i x )lj 

+[ttf - 7 x (vt x + v?y + W? x )]k|j 
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The inviscid condition is satisfied by v tl = Va, while 
Vm = V n7 specifies the impermeable condition. Re- 
flecting total energy and density as even functions, 
then a Riemann problem can be solved at the wall to 
determine the flux. This amounts to the wall being 
represented as a contact discontinuity by constraining 
the contravariant velocity to vanish. Implicitly, this 
results in 

6Q 1 = R~ l ER\ wa nSQ2 (4) 

where 


R-'ER = 


1 0 0 

0 1 - 2s* -2? t ? v 

0 —2? y ? t 1 — 27* 

0 

0 0 0 


0 0 ■ 

—2? z ? t 0 
0 

1 - 2S* 2 0 
0 1 


The block tridiagonal system may be written as 


\B[ C[ 1 


' 6Q1 ' 


' AQx ■ 

A' 3 B' 3 C' 


6Q3 

= 

AQj 


Now the change in flux across an arbitrary cell wall 
boundary is given by A 22*0/1 = A + 6Q i + A~6Q 7i or 
the sum of the changes in the flux contribution from 
the positive and negative moving waves. Substitution 
of Eq. (4) yields 


A£w* = (A + /T l i?R + A~)6Q 7 


and it can be seen that dependence upon 8Q\ has been 
eliminated. Hence, the block tridiagonal system may 
be represented with embedded boundary conditions 
as 


[ B" CJ 
A : l BJ C; 



' 6Q , ’ 


■ A Qx ' 


6Qs 

zz. 

AQs 


where B" is the appropriately modified Jacobian. 

The viscous impermeable wall imposes additional 
constraints on the specification of the wall flux. 
Again utilizing a primitive variable vector V = 
\p,u,v,w,e] T y then 

dV x 

SVi = diag[ 1, -t ut -t,, -i v , -t]SV 7 = Qy^ SV * 


for a wall face at |. In this form the toggles 
and t are set at -1 or 1 for a slip or a no- 
slip condition, or adiabatic or isothermal wall, respec- 
tively. This may be seen by simply rearranging ex- 
pressions of the form = |(ui + v a) or UvaH = 

uj u 7 . Having already specified the impermeable 
wall conditions earlier, only the viscous terms at the 
wall are of present concern. Looking at the terms of 
the form 

A Q 2 = ^{-MNiSQx + MN26Q2 + • • •) 


then substitution of the wall relations above leaves a 
term 


A <?, 


At 

V 3 


,,dVi dv 2 

Q,+ w, Q,+ 


S|«('-8) 


dv 7 

d<h 


8Q 2 + 1 


which is subsequently embedded into the block struc- 
ture. The dependent variables within cell 1 are spec- 
ified according to boundary-layer theory, holding the 
pressure gradient zero normal to the wall. The re- 
maining variables follow from fluid and thermody- 
namic relations. 

The class of problems investigated here revealed 
that the use of block implicit boundary conditions 
resulted in significantly enhanced convergence. This 
beneficial effect is caused by the faster signal propa- 
gation arising from the incorporation of the boundary 
conditions within the linear system. 


Results and Discussion 

The methods introduced in the previous sections are 
applied to test cases which demonstrate the capabili- 
ties of the algorithm. The viscous term treatment in a 
low Mach number regime is shown in the Couette flow 
problems, which are compared to similarity solutions 
and previously obtained numerical results. Demon- 
stration of the inviscid term treatment is shown by 
the capturing of transient discontinuities in the shock 
tunnel starting problem. The three-dimensional re- 
sults are compared with an experimental study of a 
hemicylinder mounted in a shock tube. 

The Couette flow problem is used to compare the 
present method against the method of Beam and 
Warming 25 and the similarity solution as given by 
Schlichting. 26 The solutions shown in Fig. 4 were 
obtained using quiescent initial conditions and vis- 
cous boundary conditions with no-slip adiabatic walls. 
Both of these cases were implemented in the thin-layer 
form at a Reynolds number of 6.4, based on the dis- 
tance between the plates, equal to 10“ 5 feet. During 
the course of these solutions, slightly more than an 
order of magnitude drop in ||p||oo p« two subitera- 
tions was observed in the (3 x 10) cell domain. The 
Courant number used for the oscillating plate calcu- 
lation was approximately 10, indicating the viability 
of these types of unsteady computations at Courant 
numbers greater than unity. The Courant number 
was computed using v = ^maxfKU, 9, W)\ + c]||s|| 
over each cell in the domain. Identical results were 
obtained using both the two- and three-dimensional 
implementations in all directional permutations. Re- 
sults reveal steeper gradients than that of the conven- 
tionally differenced scheme or the analytic solution, a 
possible consequence of the handling of the boundary 
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conditions or the viscous term treatment. In addition, 
this case was found to be insensitive to the choice of 
the higher-order flux correction terms, possibly be- 
cause of the dominance of diffusive effects. 

The third test case evaluated the inviscid term 
treatment through the simulation of the transient 
starting process of a planar shock tunnel. The (300 x 
60) cell domain is shown in Fig. 5. The solution of 
the Euler equations is presented in Fig. 6 as a com- 
parison of the experimental and numerical shadow- 
graph images, the former due to Amann, 37 while the 
graphical presentation of the latter is due to Buning 38 
This solution was obtained using Roe flux difference 
splitting with ^ = 1/3, the upwind biased flux evalu- 
ation. The Steger- Warming flux evaluation with Roe 
averaging was found to be moderately less stable, but 
no significant differences in the results were found for 
this case. The maximum compression parameter was 
used, and the entropy-fix parameter used in Harten’s 
formulation was set to 0.15. Disconcertingly nonphys- 
ical solutions were produced for smaller entropy-fix 
levels, possibly associated with an entropy-violating 
condition. The problem was initialised with a mov- 
ing shock propagating to the right at a Mach number 
of 2.97, while the boundary conditions were specified 
as impermeable inviscid along the walls and fixed for 
the inlet and exit. For the maximum Courant number 
of four used here, four subiterations were chosen per 
time step based on a subjective judgment of discon- 
tinuity sharpness. The ||p||oo was observed to drop 
approximately an order of magnitude over the course 
of these four subiterations. 

Physically, this nozzle starting process generates a 
reservoir of more than 50 times the initial pressure, 
while the density increases by 11 times the initial 
state. This reservoir provides the energy necessary 
to generate high Mach flows in the diverging nozzle 
region for short durations. The ensuing reflections of 
the shock with the nozzle wall reveals the complexi- 
ties of the shock-shock and shock-contact interaction. 
In particular, it can be seen that the development of 
the rearward-facing shock, which is directed upstream 
while being swept downstream, is resolved. At later 
times, the finer scale fluid motion between the pri- 
mary and rearward facing shocks is for the most part 
lost because of grid coarseness and attendant numer- 
ical dissipation. However, increasingly fine structures 
are captured as the grid is refined. 

The viability of the technique in three dimensions 
is shown by the final test cases. These results are in- 
tended to replicate the conditions in an experimental 
study of a hemicylinder in a shock tube by Kingery 
and Bulmash. 39 The experiment test configuration 
and pressure transducer locations are shown in Figs. 7 
and 8. In order to estimate the costs and benefits 
of inviscid versus viscous simulations, the flow about 
this geometry was computed using both the Euler 


and Navier-Stokes equations. However, the expense 
of these three-dimensional simulations permitted the 
use of only one of the inviscid flux evaluation methods; 
the Roe flux difference splitting was chosen. 

The simulation was initialized as a planar shock 
translating at a Mach number of 1.518 before diffrac- 
tion over the cylinder began. The Reynolds num- 
ber per meter was 23.3 x 10°, computed using Too = 
288.17 K t poo = 101 325 Nfm? y and specifying the ref- 
erence velocity based on the quiescent state. Bound- 
ary conditions are specified for the viscous, single zone 
computation as follows. In the {-direction, extrapo- 
lation is used. This nonphysical extrapolation is ade- 
quate for the duration of the early interaction. How- 
ever, solutions at larger times are suspect, where times 
after the shocks have propagated through the bound- 
aries are defined as large. Additionally, the use of an 
advancing front boundary is enabled because of a pri- 
ori knowledge of the grid structure and the primary 
shock speed. This simple time-dependent boundary 
reduces the computational time requirements. In the 
17-direction, the lower boundary defines the surface ge- 
ometry of the hemicylinder, and hence is specified as 
a no-slip isothermal wall. The top of the domain in 
the 17-direction, corresponding to the inner radius of 
the shock tube, is specified as an inviscid wall, based 
on the assumption that the viscous effects on this sur- 
face have negligible influence on the results. Finally, 
the {-direction boundaries are treated using the vis- 
cous condition along the floor of the tube. Symmetry 
conditions are used along the plane running along the 
longitudinal axis of the cylinder and normal to the 
floor. To simulate experimental conditions, the wall 
temperature was set equal to the temperature of the 
quiescent flow prior to primary shock arrival. The vis- 
cous grid has normal spacing of approximately 10~ K 
meters at the viscous walls. The Euler computation 
used the inviscid boundary conditions previously dis- 
cussed where appropriate. For this Euler grid, since 
the areas of the faces corresponding to the geomet- 
ric axis singularity are zero then F • s is also zero. 
Compensation for the round-off error inherent in the 
grid was implemented by eliminating those face areas 
which fell below a specified tolerance. The grids and 
boundary conditions for these cases are partly shown 
in Figs. 9 and 10. 

The inviscid computation used <f> = 1/3, /? = 4, 
Harten’s entropy fix parameter of 10" 4 , Roe flux dif- 
ference splitting, one subiteration per second-order ac- 
curate time step, and a Courant number of 15. The 
solution was obtained in 1500 time steps without any 
change of parameters. 

The viscous computation used the same flux evalu- 
ation as above with the addition of the second-order 
accurate full viscous terms. Because of the viscous 
spacing, the Courant number utilized was 10 4 , allow- 
ing the solution to be obtained in 6800 time steps 
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with no subiterations. In contrast to the inviscid sim- 
ulation, the advancing front boundary condition was 
utilized in this case. 

Results, given in Figs. 11 through 16, show that the 
primary shock is captured over two to three cells, the 
large physical thickness obtained is an artifact of the 
coarse grid used. Adaptive gridding methods would 
help maintain sharp shocks, but the anticipated ex- 
pense of these methods precluded their use here. Fig- 
ures 11 and 12 show comparisons of the numerical 
and experimental pressure histories. It is seen that 
the peak overpressure is underpredicted by 10%, pos- 
sibly owing to the coarseness of the grid which in turn 
thickens the shock. This thickening causes the poten- 
tial energy in the form of pressure to be transferred 
to the surface over several time increments. These 
computed surface pressures were extracted from the 
domain through the use of a Newton search in three- 
space for the cell in which the given (x,y, z) probe 
coordinate fell. 38 Subsequent trilinear interpolation 
over the cell, where the uniform parametric coordi- 
nates (u, u,ti>) are determined from the positions of 
the vertices of the hex&hedral cell, allows the pressure 
to be computed. Inherent in this first-order approx- 
imation lies the assumption that over a discrete cell 
the variation of pressure is linear in space. 

Figures 13 through 16 show portions of the viscous 
simulation at selected times. Physically, the inter- 
action process begins with the normal impact of the 
incident shock with the front face of the hemicylinder. 
At this time, peak overpressures of six times, and den- 
sities of four times that of the quiescent state are gen- 
erated along this forward face. As the shock diffracts 
over the sharp corner of the target, a separation bub- 
ble forms, which eventually envelops a large portion 
of the circumferential face of the body. This vortical 
motion is depicted in Fig. 16 by instantaneous stream- 
lines. A supersonic pocket is generated as the air ne- 
gotiates the sharp corner as it rushes from the stag- 
nation region left in the wake of the upstream prop- 
agating reflected shock. The next significant event 
occurs as the shock diffracts over the rearward face, 
shedding a strong vortex sheet while an expansion 
wave propagates away in a pattern which grows with 
time. The diffracted shock then impacts the floor of 
the shock tube, reflecting it upwards, while the shock 
which diffracted over the circumferential face reflects 
inwards from the outer walls of the tube. The subse- 
quent diffractions and reflections result in the inter- 
action of shocks, expansion fans, vortices, and devel- 
oping boundary layers. From experimental evidence, 
this gross unsteadiness does not dissipate for more 
than 15 milliseconds after the interaction event begins. 
However, the primary shock passes from the test sec- 
tion 5 milliseconds after the initial target interaction; 
therefore, the computation is stopped at that time. 

The effects of the viscous terms are seen by com- 


paring the pressure histories in Figs. 11 and 12. While 
the pressures along the upstream face are largely un- 
changed, the circumferential and downstream faces 
are significantly affected by viscosity. The large sepa- 
ration along these faces causes low pressure regions 
due to this vortical motion. This phenomenon is 
more accurately captured in the viscous simulation, 
as may be seen by inspection of the pressure histories 
at probe 11. Differences between the experimental 
and the present results may be due to poor capturing 
of the vortex strength owing to grid coarseness. How- 
ever, the higher-order behaviour of the method used 
here attempts to reduce the need for finely spaced 
meshes. In addition, the occurrence of deformation 
of the shroud wall is thought to be a possible event 
during the experiment. 80 

A limited cost /benefit study of the Euler versus 
Navier-Stokes equations was also performed for the 
hemicylinder case. For approximately 5.5 ms of flow 
history on a (78 x 50 x 25) cell grid, the Euler computa- 
tion consumed 7.6 processor hours, while the viscous 
simulation required 18.2 hours. From these results, 
the somewhat more accurate solution given by the 
Navier-Stokes simulation may be worthwhile. This 
is particularly true if the flowfield behaviour after pri- 
mary shock passage is important. 

Conclusions 

The application of two upwind schemes to un- 
steady, multidimensional problems within a struc- 
tured finite- volume framework has been demonstrated 
on the viscous three-dimensional blast- wave problem. 
The use of time-conservative differencing and an ap- 
proximate Riemann solver coupled with total varia- 
tion diminishing methods has resulted in time accu- 
rate nonoscillatory flowfield resolution. Newton subit- 
erations are utilized to reduce the numerical approx- 
imations made, such as factorization error and the 
inclusion of only the first-order terms in the forma- 
tion of the inviscid Jacobian. In addition, analysis 
and application of two flux evaluation methods found 
their differences to be small. Finally, for the blast- 
wave/target interaction problem the effect of viscosity 
was increasingly significant at later times. 

The expense of these algorithms is relatively high: 
86jjs per cell per iteration using a single processor 
on the Ames Research Center CCF Cray Y-MP/832 
for these vectorized codes which have computation 
rates of approximately 142 MFLOPS. In addition, the 
memory requirement is 38 words per cell. Decreased 
processor times may be achieved by many schema. 
Freezing the Jacobian for several subiterations will 
offer a processing time reduction of 15% per subit- 
eration, albeit at the expense of memory. 

While the additional computational cost of these 
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algorithms may be justified only for cases requiring 
their capabilities, the expenditure of human effort is 
likely to be less than with conventionally differenced 
techniques because of the natural scaling of the dissi- 
pation with the eigensystem. This feature provides 
somewhat greater freedom from the numerics, per- 
mitting the scientist to concentrate on the physical 
problem and the engineer on design. 

Further efforts to increase the accuracy and effi- 
ciency of these methods may be directed along the 
use of nonfactored schemes or implementation on par- 
allel machines. Geometries of realistic complexity 
will require a sonal approach, necessarily conservative 
because of the strongly unsteady compressible flow 
regimes considered here. Efficient adaptive grid tech- 
niques will reduce the memory and time expense. Syn- 
thesis of dynamical, structural, and fluid flow effects 
may provide the capability for an interdisciplinary 
simulation of the physical processes involved in this 
class of problems. 

The development of this type of tool partially ful- 
fills the objective of augmentation of experimental test 
programs, possibly eliminating the need to test cer- 
tain specific configurations altogether. Although the 
use of the inviscid flowfield time history is feasible in 
the design process, the goal of nesting full unsteady 
Navier-Stokes methods within this cycle awaits the 
advent of computers of increased power. 
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Fig. 1 Discrete hexahedral cell and stencil. 
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Fig. 6 Shock tunnel case: Shadowgraph comparison for a planar shock with initial speed of Mach 3. 
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Fig. 7 Hemicylindei case: Experiment configuration 
(after Kingery and Bulmash). 
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Fig. 9 Hemicylindei case: Central region of the 
78 x 50 x 25 cell (a)inviscid and (b)riscid grids. 
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Fig. 8 Hemicylinder case: Pressure 


transducer locations. Fig. 10 Henucylinder case: Symmetry plane of the 
viscid grid. 
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Fig. 16 Hemicylinder case: Instantaneous streamlines 
at (a)t=1.7ms, and (b)t=6.2ms. 


Fig. 16 Hemicylinder case: Density contours 
at (a)t=l.lms, (b)t=1.7ms, and (c)t=6.2ms. 
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